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_1j' Abstract 

Background: Stochastic modeling and simulation provide powerful predictive methods for the intrinsic un- 
derstanding of fundamental mechanisms in complex biochemical networks. Typically such mathematical models 
^ ' involve networks of coupled jump stochastic processes with a large number of parameters that need to be suitably 

calibrated against experimental data. In this direction, the sensitivity analysis of the reaction network's parameters 
^^ , is a crucial mathematical and computational tool since it allows to infer information about the robustness and 

O^ ■ the identifiability of the model's parameters. However, existing sensitivity analysis approaches such as variants 

^^ ' of the finite difference method can have an overwhelming computational cost in models with a high-dimensional 

T^lj- I parameter space. 

^D ' Results: We develop a sensitivity analysis methodology suitable for complex stochastic reaction networks with 

m ' a high number of parameters. The proposed approach is based on Information Theory and relies on the quantifi- 

cation of information loss due to parameter perturbations between time-series distributions, hence referred to as 
"pathwise". This is achieved by employing the rigorously-derived Relative Entropy Rate (RER), which is directly 
computable from the propensity functions. A key aspect of the method is that an associated pathwise Fisher 
K^ ' Information Matrix (FIM) is defined, which in turn constitutes a gradient-free approach to quantify the parameter 

j_j I sensitivities. The study of the structure of the FIM which turns out to be block-diagonal reveals hidden parameter 

C^ . dependencies and sensitivities in reaction networks. 

Conclusions: As a gradient-free method, the proposed sensitivity analysis provides a significant advantage 
when dealing with complex stochastic systems with a large number of parameters. In addition, the knowledge 
of the structure of the FIM can allow to efficiently address questions on parameter identifiability, estimation and 
robustness. The proposed method is tested and validated on two biochemical systems, namely: (a) the p53 
reaction network where quasi-steady stochastic oscillations of the concentrations are observed, and for which 
continuum approximations (e.g. mean field, stochastic Langevin, etc.) break down due to persistent oscillations 
between high and low populations, and (b) an Epidermal Growth Factor Receptor (EGFR) model which is an 
example of a high-dimensional stochastic reaction network with more than 200 reactions and a corresponding 
number of parameters. 
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1 Background 

The need of an intrinsic understanding of the interplay between complexity and robustness of biological 
processes and their corresponding design principles is well-documented, see for instance [THS]. The concept 
of robustness can be described as "a property that allows a system to maintain its functions against internal 
and external perturbations" [3]. When referring to mathematical models of complex biological processes, 
one of the mathematical tools to describe the robustness of a system to perturbations is sensitivity analysis 
which attempts to determine which parameter directions (or their combinations) are the most/least sensitive 
to perturbations and uncertainty or to errors resulting from experimental parameter estimation. Recently 
there has been significant progress in developing sensitivity analysis tools for low-dimensional stochastic 
processes, modeling well-mixed chemical reactions and biological networks. Some of the mathematical tools 
included log-likelihood methods and Girsanov transformations |MH]j polynomial chaos 'S , finite difference 
methods and their variants [lOilllj and pathwise sensitivity methods [12]. However, existing sensitivity 
analysis approaches can have an overwhelming computational cost, either due to high variance in the gradient 
estimators, or in models with a high-dimensional parameter space. Such issues and comparisons between 
methods where discussed in [13 . 

The aforementioned methods focus on the sensitivity of stochastic trajectories and corresponding aver- 
ages. However, it is often the case that we are interested in the sensitivity of probability density functions 
(PDF), which in nonlinear and/or discrete systems are typically non-Gaussian. In that latter direction, 
there is a broad recent literature relying on information theory tools, and where sensitivity is estimated 
by using the Relative Entropy and the Fisher Information Matrix between PDFs, p roviding a quantifica- 
tion of information loss along different parameter perturbations. We refer to [UMTS] for the case when the 
parametric PDF is explicitly known. For instance, in |16j the parametric PDF's structure is known as it is 
obtained through an entropy maximization subject to constraints. Knowing the form of the PDF allows to 
carry out calculations such as estimating the relative entropy and identifying the most sensitive parameter 
combinations. Furthermore, the pathwise PDFs are also known in reaction networks when a Linear Noise 
Approximation (LNA) is employed and in which case the relative entropy can be explicitly computed allow- 
ing thus to carry out a sensitivity analysis, |18j . However, for complex stochastic dynamics of large reaction 
networks, spatial Kinetic Monte Carlo algorithms and molecular dynamics, such explicit formulas for the 
PDFs are in general not explicitly available. 

In [19! we address such challenges, by introducing a new methodology for complex stochastic dynam- 
ics based on the Relative Entropy Rate (RER) which provides a measure of the sensitivity of the entire 
time-series distribution. Typically the space of all such time-series is referred in probability theory as the 
"path-space" . RER measures the loss of information per unit time in path space after an arbitrary pertur- 
bation of parameter combinations. RER and the corresponding Fisher Information Matrix (FIM) become 
computationally feasible as they admi t ex plicit formulas which depend only on the propensity functions (see 
(|n]) and (|S]), respectively). In fact, in [T5] we showed that the proposed approach to sensitivity analysis has 
the following features: First, it is rigorously valid for the sensitivity of long-time, stationary dynamics in 
path space, including for example bistable, periodic and pulse-like dynamics. Second, it is a gradient-free 
sensitivity analysis method suitable for high-dimensional parameter spaces as the ones typically arising in 
complex biochemical networks. Third, the RER method does not require the explicit knowledge of the equi- 
librium PDFs, relying only on information for local dynamics making it suitable for non-equilibrium steady 
state systems. In |19j we demonstrated these features by focusing on two classes of problems: Langevin 
particle systems with either reversible (gradient) or non-reversible (non-gradient) forcing, highlighting the 
ability of the method to carry out sensitivity analysis in non-equilibrium systems; and spatially extended 
Kinetic Monte Carlo models, showing that the method can handle high-dimensional problems. 

In this paper, we extend and apply the pathwise sensitivity analysis method in |19) to biochemical 
reaction networks, and demonstrate the intrinsic sensitivity structure of the network. Such systems are 
typically modeled as population stochastic jump Markov proces ses and they are simulated using either exact 
algorithms such as the Stochastic Simulation Algorithm (SSA), [20f[22] and the next- reaction method [ 23] o r 
by employing approximations such as mean field ODEs, tau-leap |24| and stochastic Langevin methods |25j. 
We show that the proposed pathwise method allows us to discover the intrinsic sensitivities of the reaction 
network by decomposing the FIM into diagonal blocks. The block-diagonal structure of the proposed FIM 
reveals, in a straightforward way, the sensitivity interdependencies between the system parameters. For 
instance, if each propensity functio n d epends only on one parameter -usually the reaction constant- then 
the FIM is a diagonal matrix (see ([15])). The sparse representation of the FIM can be essential in optimal 
experimental design as well as in parameter identifiability and robustness where each subset of the par ame ters 
defined by a block of the FIM can be treated separately. Moreover, our earlier rigorous analysis [19] for 
the stationary process regime suggests suitable extensions in the transient case which are here tested and 
validated. Finally, we present strategies for efficiently and reliably implementing the proposed method for 
high-dimensional, complex stochastic systems using an array of existing accelerated versions of the SSA 
algorithm such as mean field, stochastic Langevin, r-leap approximations and their variants, ;21.i24,i27). 



62 We test the proposed set of methods and computational strategies in two examples of complex biochemical 

63 networks. First, we study the parameter sensitivities of a p53 gene model for cell cycle regulation and 

64 response to DNA damage, that incorporates the feedback between the tumor suppressor p53 and the oncogene 

65 Mdni2 [28 . This is a reaction network that exhibits random oscillations in its steady state, and for which 

66 continuum approximations of the SSA such as LNA break down due to persistent oscillations between high 

67 and low populations. Using the proposed method, we also study a far more complex network, the epidermal 

68 growth factor receptor (EGFR) model, describing signaling processes between mammalian cells [29-31 . 

69 This is a high-dimensional system both in the number of variables and parameters, including 94 species 

70 and 207 reactions. Having a gradient-free method such as FIM for this example with parameter space of 

71 dimension 207 provides a significant advantage over gradient methods such as finite differencing, where the 

72 computation of a very high number of partial derivatives and /or directional derivatives is needed and with 

73 possibly significant variance that scales with the dimension, [llj . By contrast, the eigenvalue/eigenvector 

74 analysis of the proposed FIM identifies the order from least to most sensitive directions (determined by the 

75 eigenvectors of the FIM) by correspondingly ordering the eigenvalues. 

76 In Methods we present the derivation of the Relative Entropy Rate and its corresponding Fisher Informa- 

77 tion Matrix for continuous-time jump Markov processes as well as we reveal the block-diagonal structure of 

78 the FIM for commonly encountered reaction networks, continued by the presentation of both unbiased and 

79 biased -but accelerated- statistical estimators for RER and FIM. Then in the Results we applied and vali- 

80 dated the proposed pathwise sensitivity analysis methodology in two complex biological reaction networks. 

81 Methods 

82 We consider a well-mixed reaction network with N species, S = {^i, ..., S'jv}, and M reactions, R = 

83 {Ri, ...,i?Af}- The state of the system at any time i > is denoted by an A^-dimensional vector X(i) — 

84 [Xi{t), ...,Xjv(t)]"^ where Xi{t) is the number of molecules of species Si at time t. Let the A^-dimensional 

85 vector Vj correspond to the stoichiometry vector of j'-th reaction such that Vi_j is the stoichiometric coefficient 

86 of species Si in reaction Rj. Given that the reaction network at time t is in state X(i) == x, the propensity 

87 function, aj(x), is defined such that the infinitesimal quantity aj{'x.)dt gives the transition probability of the 

88 j'-th reaction to occur in the time interval [t, t + dt\. Propensities are typically dependent on the state of the 

89 system and the reaction conditions (i.e., external parameters) of the network such as temperature, pressure, 

90 etc. 

91 Mathematically, {X(i)}igR^ is a continuous-time, time- homogeneous, pure jump Markov process with 

92 countable state space E C N^. The transition rates of the Markov process which completely define the 

93 process are given by the propensity functions aj(-), j — 1, ..., M. Indeed, the transition rates determine the 

94 clock of the updates (or jumps) from a current state x to a new (random) state x' through the total rate 

95 ao(x) := X^fci '^j(^) which is the intensity of the exponential waiting (or sojourn) time for a jump from 

96 state x. The transition probabilities for the embedded Markov chain {Jnj ^„ with J„ := X(i„) provided 

97 the current state x are pj (x) — °'^y^' , There are exact algorithms for the simulation of the reaction network 

98 such as stochastic simulation algorithm (SSA) of Gillespie [201l32| or next reaction algorithm of Gibson and 

99 Bruck '23i as well as approximation algorithms such as r-leap [24] and several variations of it [26j[27]. As a 

100 demonstration, given that the system is at the state X(i) = x at time t, SSA computes the waiting time 5t 

101 as a random number drawn from an exponential distribution with the total rate ao(x) as parameter while 

102 the Rj* reaction occurs where j* £ {1, ..., M} is chosen such that J2'j=i Pji^) < u < X^-i^i* Pji^) where u is 

103 a random number uniformly chosen in the interval [0, 1]. The new state is given by 'K.(t + 5t) = x' = ji + Vj* . 

104 Relative Entropy 

105 Assume that two probability distributions (or more generally probability measures) V and V have cor- 

106 responding probability densities p — p{x) and p = p{x). Then, the Relative Entropy or Kullback-Leibler 

107 divergence of V with respect to V is defined as [33ll34] 

^(^i^)-A(-)i°s(iy)^- (1) 



108 In a more general setting, relative entropy is defined a,s TZlVlV) := / log ( ^ j dV where ^ is a function 

109 known as Radon-Nikodym derivative while the integration is performed with respect to the probability 
no ■p, [35 . Relative entropy has been utilized in a diverse range of scientific fields from statistical mechanics 

111 to coding in telecommu nica t ions (information theory) and finance, and it possesses the following three 

112 fundamental properties, [M], [55] : 

113 (i) it is always non-negative, 

114 (ii) it is equal to zero if and only ii V = V P-almost everywhere, and, 

115 (iii) <TZlV\'P) <ooif and only if V and V are absolutely continuous with respect to each other. 
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From an information theory perspective, relative entropy measures the loss of information when V is utilized 
instead of V, [M]. In other words, relative entropy measures the inefficiency of assuming a perturbed 

distribution V instead of assuming the true distribution V. Therefore, even though not a metric, relative 
entropy is a suitable quantity for the assessment of parametric sensitivity analysis since the higher the relative 
entropy due to some perturbation the larger the information loss in this direction. 

Moreover, relative entropy can provide an upper bound for a large family of observable functions through 
the Pinsker (or Csiszar-KuUback-Pinsker) inequality. The Pinsker inequality states that the total variation 

norm between V and V is bounded from above by the square root of the twice of the relative entropy |34j . 
However, we are usually interesting on observable functions such as mean populations, variances, correlations, 
etc. where the Pinsker inequality states that for any bounded functions / holds: 



J2n(v\vy 



\Ev[f]~E^[f]\<\\f\\oo^2n[P\V). (2) 

An obvious conclusion that is immediately drawn from the above inequality is that if the (pseudo-)distance 
between two distributions as measured hy TZ IV \V] is controlled, then the error between any two observables 
is also controlled accordingly. 

Path wise Relative Entropy and Relative Entropy Rate 

Proceeding to the pathwise formulation of the relative entropy, we assume that the propensities functions 
depend on a parameter vector 9 £ R^ (i.e., aj(x) = a^(x)) while the Markov process {X(i)} lies in 

the stationary regime. We denote by /i^(x) the steady state (or stationary) distribution of the stochastic 
process X(t). The stationary path distribution of the process in the interval [0,T] is denoted by QfoTl- 

Notice that path distributions (i.e., time-series distributions) are high-dimensional complex objects; for 
instance, if we consider the simpler discrete-time Markov chain case {X„}„gz+, defined by the transition 
probability density j)(a;, x'), then, utilizing repeatedly the Markov property, the stationary path distribution 
of a time-series (path) (xq, a;i, . . . , xt) is given by 



Q[o,T]({X„ =x„}o<„<t) =Prob(xo,...,XT) = /i(xo)p(xo,Xi) . . .p(xt-i,xt) . (3) 

138 Furthermore, we consider the jump Markov process {X(t)}fgR^ defined by perturbing the propensity func- 

139 tions by a small vector e G M.^ . The corresponding steady state and path distributions of {X(i)}4gH_^ are 

140 denoted by /x^"'"'^(x) and Qffl'^i, respectively. Let the two path distributions Qf^ j,, and Qfo^^i be absolutely 

141 continuous with respect to each other which is satisfied when a^(x) = if and only if aj~^'^{x.) — holds for 

142 all X e _E and j = 1, ..., M. Then, the Relative Entropy of path distribution Q?q j,-, with respect to Qrn j'l is 



defined similarly to ([TJ as 



^ ( Q[0,T] I Qto%) ■■= /l«g {^^] ^^fo.T] , (4) 



144 where L"^^^ is the Radon-Nikodym derivative of Qfg j,, with respect to Qfot^i • In fact, by Girsanov's formula, 

145 we can obtain an expHcit expression for the Radon-Nikodym derivative in terms of the propensity rates, [35]. 

146 In the context of sensitivity analysis, the higher the relative entropy towards a direction e, the larger the 

147 information loss t this direction, thus, the path distribution is more sensitive at this direction. Additionally, 

148 time-series distributions contain all the information about the stochastic process including the steady state 

149 distribution, as well as any type of observable, thus, pathwise relative entropy can be viewed as a "conser- 

150 vative" bound for sensitivity analysis. Indeed, having also in mind the Pinsker inequality ([2]), if relative 

151 entropy is insensitive to a parameter direction then any observable / is also insensitive to this direction. 

152 Moreover, in the stationary regime, relative entropy increases linearly with time, hence, Relative Entropy 

153 Rate (RER) 

n {Q' I Q'+^) := ^lim ^TZ (q%^^^ \ Qf+^j) (5) 

which is the time average of the relative entropy is well-defined, [36] . Due to (O, H (Q^ | Q^^*^) is an 
appropriate (time-independent) quantity for the measurement of sensitivity: it measures the rate of the loss 
of information due to an e-perturbation of the model parameters, in the long-time, stationary dynamics 
regime of the stochastic process, as first proposed in '19 . Furthermore, RER admits an explicit formula 
given by 

-A/ S ( \ 

n {Q' I Q^+^) ^ E^. [E 4 W l^g ^T^ - («o(x) - ar^(x))l . (6) 

Thus, from a practical point of view, RER is an observable of the stochastic process which can be computed 
numerically as an ergodic average in a straightforward way since it only requires the known propensity func- 
tions. Nevertheless, in order to carry out the sensitivity analysis at the parameter vector ^, the computation 
of RER for different e's is necessary, which for a high-dimensional parameter space may be computationally 
challenging. Thus, a sensitivity analysis methodology which does not depend on e's -such methods are called 
gradient-free- is desirable and is developed next. 
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Pathwise Fisher Information Matrix 

Even though not directly evident from ((6)), RER is locally a quadratic function of the parameter vector 
6 G M^. Indeed, RER is non- negative when e 7^ and equal to zero when 6 = thus the linear order term 
is zero. Therefore, RER is wri tten under a smoothness assumption on the propensity functions with respect 
to the parameter vector 9 as, [19]: 

H {Q' I Q'^+^) - ie^F«(g^)e + 0{\e\^) (7) 

where F-h(Q^) is a. K x K matrix that can be considered as a pathwise analogue for the classical Fisher 
Information Matrix (FIM). Indeed, as in the classical FIM for parametrized distributions [53], ([7]) is the 
Hessian of the RER which geometrically corresponds to the curvature around the minimum value of the 
relative entropy rate. The pathwise FIM contains up to second-order accuracy all the sensitivity information 
for the path distribution at 9 for any perturbation direction e. This practically means that the computation 
of FIM is sufficient up to second order for the evaluation of all the sensitivities of the path distribution 
around the parameter vector 9. 

The FIM can be utilized not only for the estimation/approximation of RER via ([7]) but also to infer 
intrinsic knowledge of the system's sensitivities |16,37;. The spectral analysis of the FIM at 9 reveals the 
(local) most/least sensitive directions of the system. Indeed, by ordering the eigenvalues of the FIM as 

it can be inferred that the most sensitive direction corresponds to the eigenvector with eigenvalue A^ while 
the least sensitive direction corresponds to the eigenvector with eigenvalue Af^ . Additionally, the classical 
FIM is one of the most useful tools for optimal experimental design. Many of the optimality criteria such as 
D-optimality where the determinant of the FIM is maximized or A-optimality where the trace of the inverse 
of the FIM are based on FIM, [37 . Even though the proposed pathwise FIM is not derived from a statistical 
description of measured data, meaning that the measurement noise is not taken into account, but from 



the stochastic model itself, it can be utilized for the optimal design of experiments to account the intrinsic 
properties of the system. In the same direction, robustness of the system to parameter perturbations or 
errors as well as parameter identifiability can be studied utilizing spectral analysis of the proposed FIM. 
For instance, parameter identifiability is satisfied when all the eigenvalues of the FIM are above a given 
threshold, 18 . 

Moreover, we have an explicit formula for the pathwise FIM which is given by 



F„(g^):=E, 



M 

E 



a^"(x)Ve logaf (x)Ve logaf (x)^ 



(8) 



The implications of this explicit formula are twofold. First, it reveals that for many typical reaction networks 
the FIM has a special block-diagonal structure which reflects the parameter dependencies of the parameters 
and it is discussed in detail below. Second, the FIM is based on the propensity functions as well as on their 
derivatives which are assumed known -actually, they define the process- thus the FIM, similarly to RER, is 
numerically computable as an observable of the process. Subsequent sections present various strategies to 
numerically estimate in an efficient fashion both the RER and the FIM. 



Sensitivity analysis on the logarithmic scale 

In many biochemical reaction networks, the model parameters differ by orders of magnitude and the 
only meaningful option for sensitivity analysis is to perform parameter perturbations relative to the overall 
magnitude. In our setup, this can be done by perturbing the logarithm of the model parameters instead of 
the parameters itself. Thus, utilizing the chain rule for V\ogef{0) = Vgf (9). \/ioge& = 0-^ef{9) where '.' 
means element by element multiplication (i.e., {a.h)^ — ak^ki k — 1^ ...,K), the logarithmically-scaled Fisher 
information matrix has elements: 



Fn{Q'°^')),i^ek0i(Fn{Q'))^, , k,l^l,...,K 



(9) 



where Yu{Q ) is given by ([8]). Similarly, the logarithmic perturbation for the RER is performed by utilizing 
the perturbation vector 6.e instead of e. Notice that ([7]) continues to be valid for the logarithmic scale. 
Indeed, it holds that 

1 



n 



{q' 



9|q9.(1+.) 



e^Y-H{Q'°n^ + 0{\e\') 



(10) 
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Block-diagonal structure of FIM 

In chemical reaction networks, reactions typically depend only on a small subset of the parameter vector. 
Mathematically, this is described as 

af(x)=a,(x;0fc,,...0fe^p, (11) 

where ki,...,k]^, G {1,...,-?^} while Lj <C i^ is the number of involved parameters in Rj. Using ([5]), it can 
be shown that the parametric dependences of the propensity functions are inherited to the FIM. Indeed, 
after grouping the reactions into subsets in such a way that each subset contains the minimum number of 
reactions having common parameters, the FIM -upon rearrangement of the parameter vector-, becomes a 
block-diagonal matrix. The FIM is then written as 



FniQ'°n 











A 



M 



(12) 



where Ai,...,A^^j are block matrices. The block matrices which are defined by the reaction subsets with 
the same parametric dependence are easily obtained by creating a graph whose nodes are the reactions and 
the parameters and the edges are their dependences. Then, the parameter nodes contained in a connected 
subgraph define a parameter subset which in turn corresponds to a block of the FIM. An illustration of 
this procedure is shown in Figure [1] where a reaction network with M — ^ reactions and K — 7 parameters 
is plotted. The parametric dependencies of the reactions are shown in the left side where 4 subgroups of 
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Figure 1: Left panel: The graph representation of the dependencies between the reactions (left column) and 
the model parameters (right column). The grouping of the parameters is then performed by pointing out 
the connected parts of the graph. Right panel: The corresponding block- diagonal structure of the FIM. In 
this example K = 7 while the largest dimension of the blocks is L == 3. 
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parameters are defined based on the graph connectivity. The resulting block-diagonal structure of the FIM 
is shown in the right hand side of Figure [1} 

Before proceeding with the theoretical computation of the FIM for various well-known classes of biochem- 
ical reaction networks, we list some of the implications of this simplified structure of the FIM in sensitivity 
analysis and elsewhere. 

(i) The sparsity of the FIM is proportional to the parametric decoupling between the reactions. Knowing 
a priori the zero elements of the FIM, there is no need to numerically compute them. It is clear that 
the computation cost for each sample drops from 0{K^) to 0{KL) where L is the largest dimension 
of the block matrices. 

(ii) The inverse of the FIM is also block-diagonal and each block of the inverse FIM is the inverse of the 
respective block of the forward FIM. This fact can be very useful for the estimation of a lower bound 
of the variance, i.e., Cramer- Rao bounds, [38, 39^ which is given by the diagonal elements of the inverse 
of the FIM. 



(iii) 



From an optimal experiment point of view, [37, and based on these new insights on the FIM structure, 
the various optimahty criteria are simphfied. Indeed, the determinant utilized in the D-optimality test 

is the product of the determinants of the block matrices (mathematically, det(F-j^) = ni=i '^^^(Ai)) 
while the trace of the inverse of the FIM utilized in the A-optimality test is the sum of the traces of 

the block inverses (i.e., tr(F-H^^) = Yl,i=i^^{^^^))- 



(iv) 



It can provide insights on how to increase the identifiability of the parameters. For instance, the 
identifiability of the fc-th parameter in the subsequent example will be increased if there is a way to 
increase the fc-th diagonal element of the FIM. Analogously, the robustness of the system to the fc-th 
parameter in the same example will be increased if the fc-th diagonal element of the FIM is decreased 
accordingly. 

Next we discuss two specific examples of biochemical reaction networks where the explicit calculation of 
the block-diagonal FIM is performed. 



Reactions with independent reaction constants 

An important class of well-mixed reaction networks has reactions being in the general form ^^cijAj + 

l3jBj -4- ..." where Aj and Bj are the reactant species while aj and /3j are the respective number of 
molecules needed for the occurrence of the reaction. The reaction constant, Oj, is the parameter of the j-th 



248 reaction. The propensity function for the j-th reaction is given as the product between a rate constant and 

249 a function of the current state x: 

aj(x)=%,(x), j = l,...,M. (13) 

Typically, 5j(x) — x^Xn/(a!/3!), however, it can take different forms depending on the modeling of the 
physical process. Overall^ the reaction network has K ~ M parameters, while each propensity depends only 

252 on one parameter, i.e., Lj = 1 in ([Tl]) for j = 1, ..., M. Proceeding, the (fc, l)-th element of the FIM in the 

253 logarithmic scale is given by 



250 
251 



(F„(Q'°^^)),^=0,0, 



M 

5]af( fx)de, \oga%^)de, logaf(x)^ 



(14) 



254 where /i^ is the stationary distribution of the stochastic process. However, it is easily obtained for propensity 

255 functions of type P^ that dg^ loga^(x) = -g-SkiJ) where S{-) is the Dirac function, therefore the FIM is a 

256 diagonal matrix with elements given by 



(F^(Q>o,.))^^^^| VKW], Jj^ 



(15) 



257 The diagonal elements of FIM are simply the stationary average of the respective propensity function. This 

258 rather unexpected result implies that the sensitivity of a reaction constant is proportional to the equilibrium 

259 average of the respective propensity function. However, notice that even though E^e [afc(x)] = dk^^f [^^(x)] 
the diagonal elements of the FIM are not linear functions of the corresponding reaction constants because 
the steady-state distribution fj,^ , depends also on the parameter vector 9. 

Moreover, due to the diagonal form of the FIM, it is straightforward to perform eigenvalue analysis and 
make inference about the most/least sensitive directions of the reaction network. Actually, the eigenvalues 
of the FIM are its diagonal elements while the eigenvectors are the standard basis vectors of K'^. Hence the 

265 most (resp. least) sensitive parameter is obtained from the largest (resp. smallest) diagonal element of the 

266 FIM (|l51) . Furthermore, P^ demonstrates that the (local) robustness of the reaction network to a specific 

267 parameter is inversely proportional to the average propensity of the corresponding reaction. Finally, another 

268 observation stemming from psp is that the larger the average of the propensity function the more frequent 

269 the respective reaction is fired, thus, it can be said that the faster reactions are more sensitive in a pathwise 

270 relative entropy sense. 

271 Michaelis-Menten kinetics 

Another important class of reaction networks is the Michaelis-Menten kinetics. This chemical network 
contains two reactions between species A and B (i.e., A -(^ B) with propensity functions given by 

"feW = -5 ; ai^d fl/c+ir 



272 This reaction sub-network which is derived under the quasi-steady-state assumption is one of the best- 

273 known models of enzyme kinetics in biochemistry. The parameter 9k (usually denoted by Vmax) represents 

274 the maximum rate achieved by the system, at maximum (saturating) substrate concentrations while the 

275 Michaelis constant 0/c+i (usually denoted by Km) is the substrate concentration at which the reaction rate is 

276 half the maximum value. The propensities of this Michaelis-Menten sub-network depend on two parameters 

277 (ife = L/c+i = 2 in (dJ)) thus the corresponding FIM block is a 2 x 2 matrix. The elements of the FIM 

278 matrix again for the logarithmic scale are given by 

( E^.K(x)-f4+i(x)] , l = k 

(F«(g'°^^)),, = -E,«K(x),^^ + <+,(x),;;^], l = k+l (16) 

[ , l^k,k + l 

279 for the fc-th row while the k + 1-th row is given by 

[ , l:^k + l,k . 



280 Strategies for the statistical estimation of RER and FIIVI 

281 Previous sections introduced and justified RER and FIM as appropriate observables for measuring the 

282 sensitivity analysis of ttie reaction network's parameters. Tfiis section presents various strategies how to 

283 efficiently estimate these quantities as ergodic averages of the stochastic process. 



284 
285 



Unbiased Statistical estimators 

Since the stationary distribution, /i^, is usually not known, both FIM and RER should be estimated 
numerically as ergodic averages. Indeed, the statistical ergodic estimator for RER is given by 



_. n-l M 

i=0 j=l " 






(x.) 



(ao(x.) 



,e+£ 



(x.)) 



(18) 
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where Sti is an exponential random variable with parameter the total rate, aQ(xi), while T — '^^^iSti is 
the total simulation time. The sequence {xi}2=Q is the embedded Markov chain with transition probabilities 

from state x^ to state x^+i is given by p^{xi) — l, \ . Notice that the weight 5ti at each step which is 

the waiting time at state x^ is necessary for the unbiased estimation of the observable |32) . Similarly, the 
unbiased estimator for the FIM is 



^ n-l M 

^n = y E '^*» E «'(^OVe loga^«(x,)Ve loga^^(x,)' 



(19) 



z=0 



J = l 
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Notice that the computation of the local propensity functions a^(xi) for all j = 1, ..., M is also needed for the 
simulation of the jump Markov process when Monte Carlo methods such as SSAs [32] is utilized. Thus, the 
computation of the perturbed transition rates is the only additional computational cost for the numerical 
RER while the additional cost for the estimation of the FIM is the computation of the derivatives of the 
(log-)propensities. Algorithm 1 summarizes the numerical computation of RER and FIM utilizing SSA for 
the simulation of the jump Markov process. 

Algoritlim 1: SSA-based numerical computation of RER and FIM. 

1. Initialize; x = xq, T == 0, H = and F = 0. 

2. FORi = l,...,n 



301 
302 
303 
304 



(a) Compute: {a|(x)}^!^^, a^(x). Compute also {a^*+'(x)}^" ^ (only for RER) and {Ve logaf (x)}j^£i 
(only for FIM). 

(b) St = log{ui) / al{-x) where ui - W(0, 1). 

(c) Update time: T = T + St 

(d) Update RER: % = % + StWf^^ a^«(x) log 4^^ - (a^(x) - 4+'{^)) . 



M 



(e) Update FIM: F^F + St X) li a|(x)V9 log <(x) Ve log a«(x) 



^AI 



(f) Find j* such that J2'j=i '^i(x) < '"200 (x) < J2jLj* aj(x) where U2 ^ U{Q, 1) 

(g) Update state: x = x + i^j. . 

3. Normalize: H = U/T and F = F/T. 
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310 Accelerated statistical estimators 

311 A typical feature of biochemical systems is that the modeled reaction network is large with hundreds or 

312 thousands of reactions with different time scales stemming from the orders of magnitude difference between 

313 the reaction rates and/or between the species concentrations, making the SSA extremely slow. A large 

314 number of multi-scale approximations of the original SSA have been developed in order to handle such 

315 issues resulting to accelerated simulation algorithms. For example, mean-field approximation ignores the 

316 fluctuat ions of the stochastic process and yields a deterministic ODE system for the mean population of the 

317 species [40l|41]. Stochastic corrections to the mean-field model such as stochastic Langevin [25] and linear 

318 noise approximation |42| can be applied in order to improve the accuracy of the simulati on. An alternative 

319 approximation of the jump Markov process is the tau-leap method proposed by Gillespie 24i where a batch 

320 of events occurs at each time-increment, r. Several improvement s on th e basic tau-leap algorithm on how to 

321 select adaptively the r [43] or on avoiding negative populations [27l|44] have been proposed, however, their 

322 performance is heavily model dependent. 

323 In this subsection, we suggest applying such accelerated approximations in order to efficiently compute 

324 the FIM and/or RER observables, while maintaining controlled bias in the statistical estimators. As an 

325 illustration, we present the well-established and well-known to biologists due to its computational efficiency 

326 mean-field approximation. To proceed, the stochastic process is written as 

X(i) = x{t) + 7^^{t) (20) 



where x{t) is the deterministic part (mean) of the process, ^(t) is the stochastic zero-mean part while rj is the 
amplitude of the stochastic term. The amplitude of the stochastic term is proportional to the inverse square 
root of the reactant populations [25i42j.45^. Thus, for large populations, the fluctuations of the time-evolving 
species populations become vanishingly small compared to the deterministic contributions. Consequently, 
the dominant part of the process is the deterministic term whose dynamics are given by the ODE system 

M 

i^{t) =^Vl,^a'^Jix{t)) , ^ = l,...,iV. (21) 

The ODE system (|2ip is also known as reaction rate equations |25.j . Restricted for simplicity to the special 
case presented in Section [I] the diagonal elements of the FIM are approximated as 



1 " 

i=l 
1 " 

^7^Y.^Uai{x{t,)+'nati)) (22) 

i=i 

1 " 

-^<5i,4(x(i0) + O(ry) 



T 



334 where we made use of (|^D|) . Moreover, the ODE system is assumed to be solved with an adaptive time-step 

335 numerical integrator up to final time T = X]i=o'^^^- Thus, for large species populations (IS";! 3> 1), the 

336 following numerical estimator for the FIM's diagonal elements 



1 

1 
=1 



{^n)k,k = TT.^t^^i(^(*^)) ,k = l,...,K (23) 



is derived. Relation (|23p suggests an algorithm similar to Algorithm 1 for the numerical computation of the 
FIM adapted, of course, to the deterministic case where instead of using SSA, an ODE solver is utilized. 

Remark 1: All multi-scale approximations try to overcome the problem of slow evolution of the system 
due to large number of species as well as the stiffness of the reactions. They are valid for large populations 
and relatively simple systems which do not exhibit complex dynamics such as bistability or intermittency. 
Indeed, as large deviation arguments [AQ' or even explicitly available formulas for escape times [AT ^ show that 
stochastic approximations cannot always capture correctly exit times, rare events, strong intermittency, etc. 
344 However, in order to simulate large biochemical systems there is no other alternative than such approximate 
models, which nevertheless need to be employed with the necessary caution. 
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Remark 2: In biochemical systems, scientists are interested not only on the steady state (i.e., stationary) 
regime (equilibrium or not) but also on the transient regime. For instance, signaling phenomena which 
crucially affect and determine the long time states of the system belong to the transient regime. The extension 
of the proposed sensitivity analysis method to the transient regime is justified by the fact that the time- 
normalized relative entropy can be also decomposed as a sum of simple integrals which results to the fact that 
the statistical estimators (jlSp and (J19p remain correct and valid. Subsequent section extensively presents 
an example of a biochemical system (EGFR) which exhibits transient behavior, and where the proposed 
sensitivity analysis tools are tested and validated. The rigorous mathematical derivation of the relative 
entropy rate for the transient regime is out of the scope of this publication and a dedicated mathematical 
article on the time dependent relative entropy rate will follow. 



Results 

The p53 Gene Model 

The p53 gene plays a crucial role for effective tumor suppression in human beings as its universal inac- 
tivation in cancer cells suggests [28,48,49]. The p53 gene is activated in response to DNA damage and it 
constitutes a negative feedback loop with the oncogene protein Mdm2. Models of negative feedback are capa- 
ble of oscillatory behavior with a phase shift between the gene concentrations. Here, we perform sensitivity 
an alys is to a simplified reaction network between three species p53, Mdm2-precursor and Mdm2 introduced 
in |28) . It consists of five reactions and seven parameters provided in Table [1] The nonlinear feedback regu- 
lator of p53 through Mdm2 takes place in the second reaction while the remaining four reactions fall in the 
special class where each reaction depends on one parameter. Due to these mechanisms a nontrivial steady 
state regime exists and can be characterized by stochastic oscillatory behavior, see for instance Figure [5] 
Hence the proposed RER and FIM methodology is directly applicable, and the corresponding FIM consists 
of 5 diagonal blocks with respective size 1x1,3x3, 1x1, 1x1, 1x1. Furthermore, the sensitivity analysis 
of this model has been performed earlier in |18j based on a linear noise approximation. Here we present a 
detailed comparison between the two sensitivity analysis methodologies, since the one proposed here does 
not involve any approximation of the stochastic network dynamics. 

Table 1: The reaction table with x corresponding to p53, yo to Mdm2-precursor while y corresponds to 
Mdm2. The state of the reaction model is defined as x = [x, yo, y]'^ while the parameter vector is defined as 
9 = [bx, ax,ak,k, by, gq, ay]'^ . 



Event 


Reaction 


Rate 


Rate's derivative 


Ri 


0^x 


ai(x) = bx 


Veai(x) = [1,0, 0,0, 0,0,0]^ 


i?2 


a; -;► 


a2(x) - a:xX+ ^^|a; 


Vea2(x) = [0, X, xy/{x + fc), -akxy/{x + fc)^ 0, 0, 0]^' 


i?3 


X ^ x + yo 


a3(x) = byx 


Vea3(x) = [0,0,0,0,a;,0,0]^ 


i?4 


yo^y 


a4(x) = aoyo 


Vea4(x) = [0,0, 0,0, 0,2/0,0]' 


i?5 


y^9 


a5(x) ^ Uyy 


Vea5(x) = [0,0,0,0,0,0,y]' 



Figure [3] shows the time-series of the species for the parameter values in Table [5J Evidently, oscillatory 
behavior is observed at this parameter regime, where persistent random oscillations occur, ranging between 
high and low populations. On the other hand, the frequency of the oscillations is less variable as it has 
been already reported both experimentally and numerically |28| . Another interesting observation is that 
the concentration of p53 species usually hits the lower bound of its admissible value (populations cannot 
be negative) which results in stochastic effects far away from Gaussianity, as can be readily seen also in 
Figure [51 

iO,XTO,k, k,by,aQ,ay 
RER as well as the FIM in the logarithmic scale are computed utilizing Algorithm 1. Logarithmic sensitivity 
analysis is preferred because the range of the parameters values varies by orders of magnitude as can be 
seen in Table [2] The upper plot in Figure |3] shows the RER as a function of time for various perturbations. 
Viewing RER as an observable, it is striking the speed of relaxation of the estimator. Within two or three 
oscillation periods RER has been converged to its value even though the three species have significant 
oscillations and stochasticity, as Figure [5] shows. A major reason for the fast relaxation is the numerical 
estimator of R ER w here the summation is over all reactions even though only one reaction takes places at 
each jump (see (HH)). Having the important property of quick convergence, global sensitivity analysis, where 



In this case, we denote by = \bx, a^jak, k, by, Oq, tty]'^ the parameter vector. The numerical estimator for 
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not only a point of the parameter regime but also large subsets of the parameter space, can be efficiently 
performed [T3] . The lower panel of Figure [3] shows the RER when only one of the parameters are perturbed 
by +10% or by -10%. Additionally, the RER computed from the FIM, utihzing 0, is also provided. The 
FIM approximation of RER is a second order approximation in terms of |e|, hence the computation of FIM is 
typically enough to fully resolve the local sensitivities of a model. Evidently, the most sensitive parameters 
here are bx and a^ while the least sensitive parameters are a^ and k. 

Table 2: Parameter values for the p53 model. 



Parameter 


6. 


Qx 


Ofc 


k 


by 


ao 


ay 


Value 


90 


0.002 


1.7 


0.01 


1.1 


0.8 


0.8 




100 



Figure 2: Molecule concentration of p53, Mdm2-precursor and Mdm2. Concentration oscillations as well as 
time delays (phase shifts) between the species are present due to the negative feedback loop. Furthermore, 
the concentration of p53 periodically approaches zero and since negative concentrations are not allowed, the 
stochastic characteristics of p53 are far from Gaussian. 



Comparison with tiie LNA-based sensitivity approach 

In [18^, the authors suggested a linear noise approximation for the stochastic evolution around the non- 
linear mean field equation, and based on this approximation a system of ODEs is derived for the mean and 
the covariance matrix of the approximation process. Since the noise of LNA is Gaussian, the mean and 
the covariance matrix contains all the sufficient information about the approximate model. Then, the FIM 
is derived and based on it, the sensitivities for each parameter are computed. Although there are regimes 
where this approximation is applicable (short times, high populations, systems with a single steady state, 
etc.), for systems with nontrivial long-time dynamics, e.g. metastable, it is not correct as large deviation 
arguments |46j show, or even explicitly available formulas for escape times |47j . Similar issues with non- 
gaussianity in the long time dynamics arise in stochastic systems with strongly intermittent (pulse-like) or 
random oscillatory behavior [50) . In the p53 model considered in J18j which has the same parameter regime 
as here, Figure [2] reveals that the time-series of the p53 populations persistently fluctuates between high and 
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Figure 3: Upper panel: RER in time for the parameter perturbation of bx (blue), k (green) and Uy (red) by 
-1-10% (i.e., £0 — 0.1). The relaxation time of the RER as an observable is very fast. Lower panel: RER for 
various perturbation directions computed cither directly (blue and red bars) or based on FIM (green bars). 



406 low values, thus the LNA approximation may not be accurate at least when the concentration of the species 

407 is low. 

408 At a coarse level, where the parameters are grouped into two classes depending on their sensitivities, the 

409 two sensitivity analysis approaches produce the same results. Indeed, by visual inspecting the lower plot of 

410 Figure [3] in the current publication and the Figure 3 in |18| , the sensitive parameters in both methods are 

411 bx, by, ak, do, ay while the insensitive parameters are ax, k. However, upon closer inspection, the two methods 

412 produce very different results. Figure |4] shows the proposed FIM (left) based on the exact (without dynamics 

413 approximations) pathwise relative entropy theory, as well the FIM proposed in |18, which is derived from 

414 the LNA of the reaction system. The results are completely different and the proposed FIM based on path 

415 distributions is sparse as expected. A striking difference between the two sensitivity methods is that in our 

416 proposed method the sensitivity of parameter bx is high while in the LNA-based method it is not (compare 

417 Figure U (dark blue) in this publication as well as Figure 3 in [H])- As a validation between the methods, we 

418 perturb bx as well as by by the same amount and observe the time-series of the species. Figure [5] shows the 

419 time-series for the three species of the model for the unperturbed parameters (blue lines), the perturbation of 

420 only bx by -f-10% (red lines) as well as the perturbation of only by by -1-10% (green lines). The same sequence 

421 of random numbers was used in all the simulations thus the time-series are visually comparable. It is evident 

422 that bx is more sensitive than by as our sensitivity analysis method predicted while the LNA-based method 

423 suggested the reverse order of sensitivity. An explanation of the performance of the LNA-based sensitivity 

424 analysis stems from the fact that the p53 species does not have Gaussian noise when the population is close 

425 to zero and this can happen quite often as Figure [2] (blue line) shows. Additionally, notice that both bx and 

426 by affect the concentration of p53 through the associated reactions thus their sensitivities are heavily biased 

427 due to the wrong statistical approximation of the p53 species. 
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LNA-based FIM 





Figure 4: The proposed pathwise FIM (left) based on RER as well as the (scaled) FIM based on LNA 
computed from the StochSens package [51| . Evidently, the proposed nrcthod uncouples the parameter 
correlations since most of the off-diagonal elements are zero. 
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Figure 5: Time-series of the species in the p53 model for the unperturbed parameter regime (blue), when b^ is 
perturbed by -1-10% (red) as well as when by is perturbed by -1-10%. The same sequence of random numbers 
was used in all simulations. The visual comparison between the time-series suggests that the time-series 
properties are more sensitive to b^ than to by. 
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Epidermal Growth Factor Receptor model 

The EGFR model is a well studied system describing signaling phenomena of (mammalian) cells 
As its name suggests, EGFR regulates cell growth, survival, proliferation and differentiation and plays a 
complex and crucial role during embryonic development and in tumor progression [52l l53|. In this paper, we 
adopt the reaction network for the dynamics of EGFR developed by Schoeberl et al. [5T] which consists of 
94 species with 207 reactions. Figure [5] presents at an abstract level the EGFR reaction network. Initially, 
the extracellular binding of EGF with the EGF receptors induce receptor dimerization. Then, two principal 
pathways. She-dependent and Sch-indcpendend, arc initiated leading to activation of Ras-GTP. Subsequently 
phosphorylation of MEK kinase through the activation of Raf kinase occurs leading to the phosphorylation 
of ERK kinase which regulates several proteins and nuclear transcription factors inside the cell. The detailed 
graphical description of the reaction network can be found in the Figures 1 & 2 of supplementary information 
in |31| . For the sake of completeness of this publication, all the reactions with their rates are provided in 
the supplementary files. The propensity functions for the reactions Ri, ..., Rgr, Rioo, ...,i?207 of the EGFR 
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Figure 6: Building blocks of the EGFR reaction network. Each module communicates with the previous 
or next module through few species only. Additionally, except the first module, all the other modules are 
double, one external (i.e., outside the cell surface) and one internal. 



440 
441 



network are written in the general form (see also (|13|)) 



aj(x) 



^j^A,-^' 



/(a!/3!), j = l,...,97,100,...,207 



(24) 



442 with the exception of reaction pair i?98, Rgg where their propensity functions are governed by the Michaelis- 

443 Menten kinetics 

aj(x) = KnaxXA,/(i^m + ^A,), J = 98, 99 (25) 

where x is the current state of the reaction system while Aj corresponds to the reacting species. The 
parameter vector where sensitivity analysis will be performed is the reaction constants. 



e = [A:i,...,fc97,Vmax,-ft'm,fcl 
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whose values are also provided in the supplementary txt file. Due to the specific values of the reaction 
constants as well as the initial population of the species (see Table [3]) , the firing rate between reactions 
differs by many orders of magnitude (highly stiff network), thus, even though there are few attempts to 
simulate the stochastic system [27] , here for the purposes of RER and FIM calculations we adopt the mean- 
field deterministic approximation in the sense discussed in the accelerated estimators subsection. We solve 
the derived system of ODEs with Matlab's routine odel5s and compute the FIM at the steady state regime 
which corresponds to the time interval [500, 700]. As Figure [5] suggests, the completion of the internalization 
process needs about 500 seconds. It should be noted here that even though the simulation of the EGFR 
is performed utilizing a deter min istic approximation model, the computed pathwise FIM has been derived 
from the stochastic network, (|14p . Therefore, under the validity of the approximation assumption (j20p . the 
computed FIM measures efficiently the sensitivities of the stochastic model in a gradient-free and sparse 
way. 



Table 3: Initial population of the species for the EGFR network. 



EGF 


EGFR 


GAP 


Grb2 


Sos 


Ras-GDP 


She 


4.98el0 


5e4 


1.2e4 


5.1c4 


6.63e4 


1.14e7 


1.01g6 


Raf 


Phosphatase 1 


Phosphatase 2 


Phosphatase 3 


MEK 


ERK 


Prot 


4e4 


4e4 


4c4 


le6 


2.2e7 


2.1e7 


8.1e4 
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Table 4: Ordering of the reaction rate constants based on their sensitivity index computed at the stationary 
regime. From left to right and from up to down. 
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kl47 


K„i 


fcl21 


kl2 


ki8 


^79 


k83 


ki9 


kis 


ku9 


kl33 


^82 


fcl76 


ki 


k22 


^25 


k87 


^206 


^34 


k35 


^200 


kl29 


kas 


kse 


k207 


k202 


^75 


fcl26 


kl30 


k201 


fc43 


kij 


^203 


fcl94 


kl89 


A:42 


kl85 


^46 


kl95 


kl91 


kl90 


fcl86 


fcl27 


kl92 


k2 


kl79 


ki7 


kl6 


fcl80 


fcl68 


^172 


kl52 


fcl55 


kl50 


kl54 


kl70 


fcl73 


fcl69 


ku 


ki7i 


kl75 


km 


kl53 


fcl51 


fcll8 


"-52 


k92 


^56 


kl09 


kl08 


kii9 


^65 


ksi 


►'max 


k73 


k53 


k32 


k33 


fcl03 


fcl02 


^41 


kue 


^50 


kgo 


^54 


kii7 


^205 


kes 


ke2 


kl98 


fcsi 


k93 


kl99 


fcioi 


^30 


^31 


kiQo 


fcl93 


kl87 


^184 


fcri 


fcl88 


fcgi 


ki8i 


kl78 


^39 


kei 


fc57 


kio 


fcli4 


^106 


keo 


^24 


ki3 


^29 


^59 


^55 


k28 


kii2 


kl04 


^23 


fe04 


fcl96 


fcl07 


kss 


ku 


kii5 


k27 
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The upper plot of Figure [7] shows the diagonal elements of the FIM in descending order computed at 
the stationary (steady state) regime. The fc-th diagonal element of the FIM corresponds to RER where the 
perturbation takes place only to the fc-th parameter (see ©). Figure [7] (upper plot) in conjunction with 
Table |4] where the reaction constants starting from the most sensitive towards the least sensitive parameter 
fully describe the (local) sensitivities of the reaction network. We report our results in the format of Figure [7] 
in order to be able to accommodate the large number of parameters in the model. Moreover, since the FIM is 
diagonal -except a small 2x2 block associated with the Michaelis-Menten reactions- the diagonal elements 
correspond to the eigenvalues of the FIM. The sensitivity analysis depicted in Figure [71 demonstrates that 
most model parameters allow for a vast range of perturbations without affecting the system's dynamics, see 
for instance Figure El Furthermore, this robustness to variations in most parameters was also reported in 
the original, fully deterministic EGFR model 31 . This is a feature shared by many multi-parameter models 
in systems biology and which is known as " sloppiness" , [54 . Our methodology can easily demonstrate such 
properties in stochastic dynamics, as we can readily see in FigurejTl even if the models include a large number 
of parameters. 

The earlier discussions refer to the analysis of the EGFR model close to steady state. On the other hand, 
EGFR is a signaling model whose transient regime, in addition to the steady state, is of great interest. As 
discussed in Remark 2 in Section [I] we can justify the application of the RER and FIM sensitivity analysis 
in the transient regime. Therefore, we compute the proposed FIM at the time interval [0, 10], using (p3|). 
The lower plot of Figure [71 shows the diagonal elements of the FIM in the transient regime while keeping 
the ordering of the parameters unchanged as in the upper plot which depicts sensitivity at the steady state. 
The parameter sensitivities are completely different meaning that the sensitivities are time-dependent in 
the transient regime. For instance, the most sensitive parameters in the stationary regime correspond to 
the final products of the reaction network however in the time interval [0, 10] these species have not been 
produced yet resulting to associated insensitive reaction constants. 

Furthermore, the Pinsker inequality ([2]) implies that the insensitive parameters can be perturbed by 
orders of magnitude without affecting the species concentrations or any other observable. As an illustration 
of this fact. Figure m presents the concentrations of various critical species of the EGFR model when the 140- 
th (fces) most sensitive parameter is perturbed, see Tabled Solid blue lines correspond to the unperturbed 
parameter case while the dashed red lines correspond to the perturbed case where the perturbation is a 
multiple by a factor of ten of the corresponding reaction constant. We chose to present the total number of 
(EGF-EGFR*)2 binding species without Sch* (top, left panel) and with Sch* (top, middle panel) as well as 
Ras-GTP (top, right panel), total activated Raf or total Raf* (low, left panel), doubly phosphorated MEK 
or MEK-PP (low, middle panel) and doubly phosphorated ERK or ERK-PP. These species are important 
for the understanding of the system since the different modules of the EGFR reaction network communicate 
critically through these species (see Figure [6]). 



16 




100 120 

Parameters 



20 r 



-20 



I -60- 
-80- 
-100 



IP 



■jkuK^A-f^V 



20 



40 



60 




100 120 

Parameters 



Figure 7: Diagonal elements of the FIM computed at the steady state regime (upper plot) and at the 
transient regime (lower plot). Presumably parameter sensitivities differ by orders of magnitude. 
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Conclusions 

In this paper, we applied and extended a recently proposed parametric sensitivity analysis methodology 
to complex stochastic reaction networks. This sensitivity analysis approach is based on the quantification of 
information loss along different parameter perturbations between time-series distributions. This is achieved 
by employing the rigorously-derived Relative Entropy Rate, which is directly computable from the propen- 
sity functions only. A key aspect of the method is that we can derive rigorously a corresponding Fisher 
Information Matrix on path-space, which in turn constitutes a gradient-free approach to quantify the para- 
metric sensitivity analysis; as such it provides a significant advantage in stochastic systems with a large 
number of parameters. We further demonstrated that the structure of the FIM revealed hidden, parameter 
dependencies and sensitivities between the reactions. The block-diagonal structure of the FIM highlighted 
the sparsity of the matrix which resulted in further increasing of the efficiency of the proposed method. 
Therefore, parametric sensitivity analysis for high-dimensional stochastic reaction systems is now tractable 
since it is well-known that in high dimensional stochastic systems the sensitivity analysis techniques can in- 
volve estimators of very high variance, e.g. in finite difference methods and their recently proposed variants, 
which can present an overwhelming computational cost |llj . Additionally, we proposed the use of multiscale 
numerical approximations of stochastic reaction networks in order to derive efficient statistical estimators 
for the FIM and tested one such approximation (mean field) to a high-dimensional system. 

The proposed pathwise sensitivity analysis method is tested and validated on two biological systems: 
(a) the p53 reaction network where quasi-steady stochastic oscillations of the concentrations are observed 
and where multiscale approximations break down due to the consistent alternation between low and high 
populations, and (b) a stochastic epidermal growth factor receptor model which is an example of a high- 
dimensional reaction network with more than 200 reactions and a corresponding number of parameters. 
In the EGFR reaction network, we combined the proposed pathwise FIM which has been derived from the 
stochastic network and the mean field approximation which is used for the efficient estimation of the pathwise 
FIM. Moreover, our earlier rigorous analysis for the steady state regime [T51 suggests suitable extensions in 
the transient regime which has been tested and validated for the EGFR model. We develop the full rigorous 
theory in an upcoming publication. 
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Figure 8: Time-dependent concentration of various species of the EGFR network either for the unperturbed 
parameter vector (sohd blue hues) or for the perturbed one (dashed red hues). The 140-th most sensitive 
parameter (fees) is ten- fold increased and the species concentrations are not affected. For the least sensitive 
parameters such as fcgS: we rigorously know form Pinsker inequality ([2|) that they should not alter the 
concentration values or any other observable even when they are heavily perturbed. 



Finally, we note that the relation between pathwise RER and various observables is not straightforward. 
However, we note that the path distribution contains all information regarding the process including the 
steady state and all time-dependent observables: practically, our proposed sensitivity analysis represents 
a conservative sensitivity estimate in the sense that insensitive directions for the relative entropy on path- 
space, will yield insensitive directions for every observable. This latter statement can be seen mathematically 
through the Pinsker inequality, ([2]). Based on these observations, the proposed sensitivity analysis methods 
can be easily used in complementary fashion with existing sensitivity analysis tools, as it can be used to 
narrow down the most sensitive directions in a system. 
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